Typeset with jpsj2.cls <ver.l.2> 



Full Paper 



Replica Analysis of a Preferential Urn Model 

Jun Ohkubo 1,2 *, Muneki Yasuda 1 and Kazuyuki Tanaka 3 

1 Department of System Information Sciences, Graduate School of Information Sciences, 
Tohoku University, 6-3-09, Aramaki-Aza-Aoba, Aoba-ku, Sendai 980-8579, Japan 

2 JSPS Research Fellow 

3 Department of Applied Information Sciences, Graduate School of Information Sciences, 
Tohoku University, 6-3-09, Aramaki-Aza-Aoba, Aoba-ku, Sendai 980-8579, Japan 

We analyze a preferential urn model with randomness using the replica method. The prefer- 
ential urn model is a stochastic model based on the concept "the rich get richer." The replica 
analysis clarifies that the preferential urn model with randomness shows a fat-tailed occupation 
distribution. We discuss the intuitive picture of the urn model, a relationship between quenched 
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1. Introduction 

Urn models have been used in order to explain various 
physical phenomena such as quantum gravity and glassy 
dynamics. 1_8 - ) Because of their tractability for analytical 
calculations, the urn models have been studied exten- 
sively in statistical physics and probability theory. The 
examples of the urn models are the backgammon model 1 -' 
and the zeta urn model; 2 ) these models can be treated 
analytically. Furthermore, it has been shown that these 
urn models relate to the zero-range process which is also 
famous in the statistical physics; see a review by Evans 
and Hanney. 9 ) 

The preferential urn model 10 ) is also one of the ex- 
amples of the urn models. It has been revealed that the 
preferential urn model shows a fat-tailed occupation dis- 
tribution only in the case with randomness; the meanings 
of the "randomness" will be explained later in the present 
paper. The fat-tailed behavior is similar to the Parato's 
distribution. 11 ) The preferential urn model stems from 
research on complex networks, 10 ) and because of the sim- 
plicity, the model would have a wide range of applications 
in research on not only complex networks, but also many 
other fields such as econophysics and social sciences. In 
econophysics or social sciences, it could be reasonable to 
consider heterogeneous cases, e.g., in which each agent 
in econophysical systems has a different ability. 

In the present paper, an analysis of an urn model with 
randomness is presented. While equilibrium properties 
of urn models without randomness are easily obtained, 6 ) 
the analysis for the urn models with randomness are not 
straightforward. We analyze the preferential urn model 
with randomness by using the replica method. 12 ) Fur- 
thermore, we will give some discussions for the analytical 
treatment and an example of the applications of the pref- 
erential urn model for econophysics, though the example 
is just a trial one. 

The outline of the present paper is as follows. Section 2 



gives a brief introduction to urn models. In § 3, we intro- 
duce the preferential urn model. The general discussion 
for calculating occupation distributions is given in § 4, 
and § 5 gives the replica analysis for the preferential urn 
model. Some discussions for the replica analysis and an 
example of the applications of the model are given in § 6. 
Finally, we draw the main conclusions in § 7. 

2. Backgrounds of Urn Models 

One of the most important properties of the urn mod- 
els is that the partition function is given by a factorized 
form. 6 ) For example, we consider a system with N urns 
and M balls which are distributed among N urns. When 
we define an energy of each urn as E(rii) by using the 
number of balls in the urn, rij, the total energy of the 
whole system is given by 

N 

E(n u ...,n N ) = J2 E ( n i)- (!) 

i=l 

The partition function of the system is calculated from 

oo oc / N \ 

Z= ■■■ J2 Pm---Pn N s(j2 ni ,M) , (2) 

ni=l riN — l \i=l / 

where 



is the unnormalizcd Boltzmann weight attached to urn 
i, and (3 is an inverse temperature. The Kronecker delta 
of eq. (2) stems from the restriction of the total number 
of balls. We can obtain various physical quantities by 
using the partition function of eq. (2); for example, we 
can calculate an equilibrium occupation probability P(k) 
with which an urn has k balls in the equilibrium state. 

The above example, i.e., eq. (2), corresponds to mon- 
key class of urn models. 6 ) In the Monkey class, an urn 
is selected at random, from which any ball is drawn, and 
transfer the drawn ball to the other urn; this dynamics 
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Fig. 1. Illustration of the preferential urn model. A ball is drawn 
at random, and then the drawn ball is moved to an urn chosen 
preferentially. 



corresponds to the image of a monkey playing at ex- 
changing balls among urns. In other words, each ball is 
indistinguishable in the monkey class, so that the equi- 
librium statistics of the monkey class is considered as 
Bose- Einstein statistics (quantum). On the other hand, 
there is the other type of urn models, so-called Ehren- 
fest class, which corresponds to the Maxwell-Boltzmann 
statistics (classical). In the Ehrenfest class, balls within 
an urn are distinguishable, so that the partition function 
of the Ehrenfest class is written as 

oc oo / N \ 



n N = l 



n N < 



which involves inverse factorials l/(n,^!); it is different 
from the partition function of the monkey class of eq. (2). 

The choice of the energy of each urn is arbitrary, and 
hence one can consider a lot of models so that a selected 
model is suitable for his/her own problem. 

3. Preferential Urn Model 

We consider a system of M balls distributed among N 
urns, and denote the density of the system as p = M/N. 
The number of balls contained in urn i is denoted by rij , 



and hence 



N 



M. For the preferential urn model, 



we define the following energy of each urn: 
E(m) = -ln(n;!). 



(5) 



From the definition of the energy, it is clear that each 
urn tends to get more and more balls because an urn is 
stable when the urn obtains a large number of balls. Note 
that there is a restriction for the total number of balls, 
so that an urn cannot get an infinitely large number of 
balls. It would be clarified that the transition rate in the 
urn model with the above definition of the energy obeys 
a preferential concept, as described later, and hence the 
urn model is called the preferential urn model. Further- 
more, we here assume that the preferential urn model 
belongs to the Ehrenfest class. 

The unnormalizcd Boltzmann weight of eq. (3) is writ- 
ten by 



Prii 



(n,!)*, 



(6) 



where fli is an inverse temperature of urn i. Note that 
each urn has own inverse temperature and each inverse 



temperature could be different from each other; we de- 
note the distribution of the inverse temperatures as 0(/3). 
Though it might seem strange in terms of physics that 
each urn has a different temperature, there should be no 
restriction for a homogeneity of temperature in informa- 
tion physics, social sciences, and econophysics (Or one 
could say that the system is in nonequilibrium.) We call 
the inverse temperature "the inverse local temperature" 
because of its locality. It will be shown in § 4 and § 5 
that the locality is important for generating a fat-tailed 
occupation distribution, and, furthermore, these inverse 
local temperatures allow us to understand results of an 
equilibrium occupation distribution intuitively. 

Using the heat-bath rule, we derive the transition rate 
Wm^m+i from the state Hi to n, + 1 as follows: 3 ^ 



Pm 



(7) 



The transition rate indicates that an urn with more balls 
tends to get more balls, and this means the preferential 
concept. The dynamics of the preferential urn model is 
illustrated in Fig. 1. 

In summary, the construction of the preferential urn 
model is as follows: 

(i) Set N urns and M balls. The M balls are distributed 
among the N urns, and the initial distribution is 
arbitrary (usually distributed at random). Each in- 
verse local temperature Pi is chosen by using a dis- 
tribution 4>(f3). 

(ii) Choose a ball at random (the Ehrenfest class). 

(iii) Move the drawn ball to an urn selected by using the 
transition rate W ni -mt+i- 

(iv) Repeat the procedures (ii) and (iii) until the system 
reaches an equilibrium state. 

If all urns have the same temperatures, i.e., (3i — [3 for 
all i, equilibrium occupation distributions are easily cal- 
culated. 2,6 ' 10 ) For instances, the equilibrium occupation 
distribution in the high-temperature limit (/% — > for all 
i) is a Poisson distribution. 6,10 ) However, it is generally 
difficult to calculate equilibrium occupation distributions 
of urn models with randomness. In the following sections, 
we give an analytical treatment to calculate the equi- 
librium occupation distribution of the preferential urn 
model with randomness. 

4. Occupation Distribution 

In this section, we give a general discussion for cal- 
culating occupation distributions, and it will be revealed 
that the case with randomness is difficult to be analyzed. 
In the next section, the replica analysis is given, and an 
explicit form of the occupation distribution will be ob- 
tained in the case with randomness. 

By using the integral representation of the Kronecker 
delta, 2wiS(a,b) = §&zz a ~ b ~ 1 , the partition function of 
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eq. (4) is calculated as follows: 



oo oo 



m— riN—0 

2m J 



(8) 



Note that all values of the inverse local temperatures 
{/3i, . . . ,Pn} are clamped at certain values in cq. (8). 
Here, we define fjf 1 ^ as an equilibrium occupation prob- 
ability of urn 1 with the inverse local temperature f3\ in 
a certain configuration: 

oo oo 



ni— n N —0 
/ N \ 



ni\ un- 



(9) 



Then, fjf 1 ^ indicates the probability with which urn 1 
has k balls in a certain configuration. Next, we define Z 2 
as 



oo oo 



z 2= E ••• E fa!)*- 1 ...^!)^- 1 

712—0 njv— 



2m 



(10) 



and simplify fjf 1 ' of eq. (9) as follows: 

/W0 = (fc! )A-i|». 

In order to obtain the equilibrium occupation distri- 
bution of urn 1 with /3i in the configurational aver- 
age, P(k,/3i), we average fjf 1 ^ over the configuration 
{02,- ■■ ,Pn}: 



(11) 



P(k,0i) = (fl™) 



{02,-,0n} 



(12) 



{02,..;Pn} 



where (• • • ){p 2t ... l p N ] means that the configurational av- 
erage is taken only over {/3 2 , • • • , 0n}; the configurational 
average of a quantity A(fli, /3 2 , • • • , (3n) is defined as 

W.ft W 

= J <A(/3 2 )d/3 2 • • ■ y 0(/3 A r)d/3 A rA(/3 1 ,/3 2 ,...,/3 A r). 

(13) 

Finally, we take the configurational average of P(k, (3) 
for [3. Then, the equilibrium occupation distribution is 
calculated from 



/ 



P(k)= / d0 4>(J3)P(k,l3). 



(14) 



In eq. (12), it is difficult to take the configurational 
average of the numerator and denominator at the same 
time. Therefore, we use replica analysis for evaluate the 
equilibrium occupation distribution. 



5. Replica Analysis 

The replica method 12 -* is one of the powerful methods 
for analysing random systems. In the replica analysis, 
we prepare m replicas of the original system, and evalu- 
ate the configurational average of the partition functions, 
assuming the thermodynamic limit. After that, the limit 
m — > is taken, which has successfully clarified typical 
properties of many random systems. While mathematical 
validation of the replica analysis still remains unresolved, 
the replica analysis has been empirically recognized as a 
reliable scheme in physics; for instance, we can use the 
replica analysis in order to calculate a free energy or or- 
der parameters. In this section, we use the replica anal- 
ysis to obtain the equilibrium occupation distribution of 
the preferential urn model with randomness. The aim 
of the present paper is to propose the analytical treat- 
ment for the preferential urn model with randomness, 
and hence we here restrict ourselves only to the case in 
which the distribution of the inverse local temperatures, 
(j)(f3), follows an uniform distribution, i.e., <p(f3) = 1 for 

pe [0,1]. 

First, we assume the following expression in order to 
calculate eq. (12): 



= exp((ln/f l) \ . (15) 

{(>2,...,M \\ fe /{/3 2 ,...,/W 

The expression is obtained from the self-averaging prop- 
erty 12 ^ because of the equilibrium occupation probability 
fjf 1 ^ is intensive. Therefore, it is necessary only to cal- 
culate the configurational average of Infjf 1 ^; this term 
is calculated from 

+ Q* Z *) {^....m- ( 16 ) 

In the replica analysis, we use the following manipula- 
tions to evaluate the configurational average: 



(lnZj) 



{02,...,0 N } J™ 



(z?) 



i !{P2,:.,Pn} 



- 1 



m 



(i = l,2). (17) 

We should calculate the following quantities to evaluate 
cq. (16) using eq. (17). Therefore, we first calculate the 
following quantities: 

( Z D{p 2 ,...,p N } 



/n{^- ( " +i)G w r " iff (^.A)}' ^ 



where 



{02,-,0n} 



n fe*« ( " +i ~ fc)G (*") jv ~ i 



(19) 



„1 oo 

G(z)= / d/3 4(0) Y,{n\f~ 1 z n 

pi oo 

= / d/3 V(n!)"-V (20) 

^0 n 
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and 



H(z,0) 



n=0 



(n!) 



\V 3 - 1 z n . 



(21) 



As mentioned above, we here treat the case with uniform 
randomness, i.e., <p{[3) = 1, /3 € [0, 1], so that the second 
line of eq. (20) can be derived. Next, we assume replica 
symmetry, i.e., z a = z, and derive the replica symmet- 
ric solution. Then, by using the saddle-point method, 
eqs (18) and (19) become as follows: 

(•^™){/32,-,/3iv} 

~ exp{— mhi2: s — mlnG(2; s ) + m\nH(z s , (3\) 
-mpN In z s + mN In G(z s ) + m In C} 

~ 1 — mlnz s — m\nG(z s ) + mh\H(z Sl [3\) 

- mpN In z s +mN In G (z s ) + mln C, (22) 
where C is a constant; 

( Z 2 l ){l3 2 ,...J3 N } 

~ 1 — m(l — A:) lnz s — m\nG(z s ) 

- mpN\nz s +mN\nG(z s ) +m\nC. (23) 

In deriving the last expression of eqs. (22) and (23), we 
consider a second-order Taylor expansion because the 
limit m — > will be taken. The saddle-points z s in 



(ZD 



r} and (ZD 



{/3 2 ,...,Pn} 



are the same and cal- 



culated by using the following saddle-point equation: 

z s G'(z s ) 



G(z s ) 



where 



G(z s ) = l + z s + Y J 



n=2 



Inn! 



1-1 



and 



G'(z s ) = 1 



E 

n=2 



Inn! 



1 

1 



(24) 



(25) 



(26) 



Equations (25) and (26) are easily obtained from eq. (20). 
From the above calculations, we obtain the following ex- 
pression by rearranging eq. (16): 



= (Pi - l)lnfc! -]nH(z s ,/3 1 ) + fclnz s 



(27) 



Therefore, the equilibrium occupation distribution of urn 
1 with pi in the configurational average is obtained from 
eqs. (12), (15), and (27) as follows: 

P(k,p 1 )^(kl)^- 1 - 



H(z s ,fa) 



{kiy 



EZo(nl)^zn- (28) 

Finally, we obtain the equilibrium occupation distribu- 
tion of the whole system as 



P(k) 



d/3 



(k\) 



(29) 




Fig. 2. Equilibrium occupation distributions of the preferential 
urn model. The numerical results for the density p = 1.0 are 
shown by squares, and those for p = 4.0 by triangles. The solid 
lines correspond to eq. (29) for respective cases. In both cases, the 
number of urns N is 1000, and that of the transfer process of balls 
is 2000 A 1 ". The data are averaged over 20 different realizations. 



Note that the equilibrium occupation distribution is suit- 
ably normalized. 

A large density p leads to z s ~ 1, and hence we con- 
sider that the saddle-point is approximately z s = 1 for 
/9> 1. Substituting z s = 1 into eq. (29), we obtain the 
asymptotic form of the equilibrium occupation distribu- 
tion as follows (see Appendix A): 

Hence, the equilibrium occupation distribution follows a 
generalized power law with a squared inverse logarithmic 
correction. 

To confirm the above analytical treatments, we per- 
formed numerical experiments. Because initial distribu- 
tions are arbitrary, we put each ball to urns at uniformly 
random as initial configurations. In the numerical exper- 
iments, the number of urns N is 1000, and that of the 
transfer process of balls is 20007V. We checked that the 
20007V transfer processes are enough to give the equilib- 
rium states. Figure 2 shows the comparison between the 
equilibrium occupation distributions of numerical exper- 
iments and those of the replica analysis. For the density 
p = 1.0, we obtain the saddle-point z s — 0.660 from 
eq. (24); z s = 0.967 for p = 4.0. The numerical experi- 
ments indicate that the replica analysis is adequate. 

6. Discussions 

6. 1 Intuitive understanding of the local temperatures 

Here, we discuss an intuitive understanding of the pref- 
erential urn model. The heterogeneity of the inverse lo- 
cal temperatures indicates that each urn has different 
ability to obtain balls. The inverse local temperature al- 
lows us to obtain the following intuitive picture; urns 
with high-temperature tend to kick out balls, and those 
with low-temperature are likely to stock balls. There- 
fore, when there is no randomness, all urns compete for 
balls, so that a fat-tailed behavior does not occur. On the 
contrary, when an urn has a lower-temperature, the urn 
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tends to get a lot of balls. Hence, it is easy to imagine 
that urns with a lot of balls can emerge due to random- 
ness. 

6. 2 Restriction for the inverse local temperature (3 
We note that the value of the inverse local temperature 

0i has a restriction: the value of can not be over 1. 
For example, eq. (20) does not converge if j3 > 1. The 
restriction of (3 < 1 means that the urn with (3 > 1 is too 
cold to kick out balls in the urn, so that the model does 
not work well. Actually, we have confirmed by numerical 
experiments that almost all balls concentrate in an urn 
with (3 > 1. Therefore, the restriction of < (3 < 1 is 
necessary. 

6.3 Annealing systems 

As mentioned in § 3, it is difficult to calculate the 
configurational average of the numerator and denomina- 
tor at the same time, i.e., the configurational average of 
cq. (12). Hence, we use the replica method to calculate 
the configurational average. The difficulty stems from the 
fact that the preferential urn model has quenched disor- 
der; evaluation of a physical quantity is taken for a given 
fixed (quenched) set of f3i, so that the inverse local tem- 
peratures are time- independent. However, if we consider 
an annealing system, i.e., the inverse local temperatures 
are not quenched and changed with time enough fastly, 
we can easily take the average of eq. (12), and there is 
no need to use the replica method. 

For the annealing system, 

(I) =f^^>, (31) 

and these quantities of (Zi)^ ^ , and {Z,^,^ ^ , 
are easily obtained. Finally, we can get the equilibrium 
occupation distribution for the annealing system as the 
same form of eq. (29). It is worth pointing out that the 
case of the annealing system gives the same results as for 
the quenched system. For the annealing system, we can 
straightforwardly get the equilibrium occupation distri- 
bution without using the replica method, then one might 
consider that there is no need to use replica method to 
analyze urn models with randomness. However, when we 
have performed numerical experiments for the anneal- 
ing system in which all inverse local temperature vary 
at each time step, the fat-tailed behavior does not oc- 
cur. Hence, it is an open question why the result for the 
annealing system agree in that for the quenched system, 
and how situation corresponds to the annealing system. 
At least, we consider that the replica method is adequate 
for analyzing urn models with quenched randomness. 

6.4 Econophysical application 

The preferential urn model stems from the analysis 
of the nongrowing models for complex networks, 10 ) but 
it would be a useful model for various phenomena. For 
example, we consider an econophysical model in which 
many agents exchange their wealth. Each agent corre- 
sponds to each urn, and the transfer procedure of their 
wealth is represented by the transfer process of balls. We 



here assume the transition rate Wk^kt+i °c (h + 1)^S 
where ki is the amount of wealth of agent i, and (3i rep- 
resents an ability of agent i for obtaining wealth. The 
above dynamics indicates the following statements: 

• A wealth to be transferred is selected at random. 
It means that agents with large wealth are likely to 
loose their wealth (Agents with large wealth have a 
lot of opportunity to consume their wealth.) 

• Agents with larger wealth have the chance to get 
wealth in higher probability. 

• Each agent has a different ability for obtaining 
wealth. 

Hence, the above model seems reasonable for an econo- 
physical one. The econophysical model can be treated 
analytically, as shown in the present paper; supposing 
that the distribution of their abilities (3i is an uniform 
one, we conclude that the wealth distribution has a fat- 
tailed form of eq. (30) when total wealth in the whole 
system is enough large. Though the fat-tailed form has 
the squared inverse logarithmic correction, the cumula- 
tive distribution of eq. (30) is similar to the Parato's 
distribution. 11 ^ The preferential urn model has not been 
proposed in order to the Parato's distribution, so that 
it is unclear that the above discussion is suitable; the 
model is just a trial one. However, it seems reasonable 
to assume that a system includes some randomness; each 
agent is in various environments and has a different abil- 
ity. Furthermore, the preferential concept also seems to 
be adequate for real econophysical systems. 

7. Concluding Remarks 

In summary, we analyzed the preferential urn model 
with randomness using the replica method. While it is in 
general difficult to analyze urn models with randomness, 
the preferential urn models with uniform randomness for 
inverse local temperatures can be treated analytically by 
using the replica method. From the analysis, it was re- 
vealed that the equilibrium occupation distribution has 
a fat-tailed form when the density is enough large. The 
fat-tailed distribution does not follow a pure power-law 
distribution; the distribution has a generalized power law 
with a squared inverse logarithmic correction. 

The preferential urn model would produce various fur- 
ther studies, such as nonequilibrium properties, other 
variation of the urn models with randomness, and ap- 
plications of the urn models for various other research 
fields. It will also be important to discuss relationships 
with interacting particle systems, such as the zero-range 
process. 9 ) 

Acknowledgment 

The authors thank Naoki Masuda for his helpful com- 
ments on the manuscript. This work was supported in 
part by grant-in-aid for scientific research (No. 14084203 
and No. 17500134) from the Ministry of Education, Cul- 
ture, Sports, Science and Technology, Japan. J. O. was 
partially supported by Research Fellowship from Japan 
Society for the Promotion of Science (JSPS) for Young 
Scientists. 



6 J. Phys. Soc. Jpn. 



Full Paper 



Author Name 



Appendix: Derivation of eq. (30) 

When we approximate the saddle-point z s = 1 for p 3> 
1, the equilibrium occupation distribution of cq. (29) is 
given by 

(fciy 3 - 1 



p(k) = r 

Jo 



d/3 



(A-l) 



By using the partial integration twice, the approximate 
equilibrium occupation distribution for large k region is 
obtained as follows: 



P(k) 



1 



ln(fc!) 



— e 



-lnfc! 



1 



H'(l,l) 



+ 



(ln(fc!))2 
1 



„-ln(fc!) H'i 1 , ) 

H(1,0) 2 



0((ln(fc!)) 3 ) 
(ln(*0) 2 (Er =0 l) 2 



(ln(fc!)) 2 i^o 
lim 



(Et i) 2 

0(S 2 ) 



(ln(fc!)) 2 s^oo 0(S 2 ) 
,_o 1 



(lnfc) 



2 ' 



(A-2) 



where we use the following relation, 13 ) 

s 

lim V ln(n!) - lim 0(S 2 ) 



S^oo 



n=0 



(A3) 



and the Stirling formula, ln(fc!) ~ fclnfc — k. 
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